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ABSTRACT 


As  part  of  the  Office  of  Naval  Research  (ONR)  sponsored  Windy  Island  Soliton 
Experiment  (WISE),  two  deep  water  moorings  were  placed  in  the  northeastern  portion  of 
the  South  China  Sea  deep  basin  to  conduct  an  acoustic  propagation  study.  For 
approximately  one  year  the  source  and  receiver  transmitted  and  received  phase- 
modulated  signals  to  measure  the  multi-scale  variability  in  the  transmission  loss  induced 
by  the  ocean  mesoscale  variability  and  the  progression  of  internal  tides  and  waves.  A 
numerical  acoustic  propagation  model  based  on  Hamiltonian  ray  tracing  is  utilized  to 
replicate  the  observed  basic  arrival  structure  and  transmission  loss.  Being  able  to 
accurately  model  the  basic  arrival  structure  is  a  necessary  first  step  before  modeling  the 
observed  variability  can  be  attempted.  The  comparison  of  the  modeled  arrival  structure 
with  the  actual  data  was  utilized  to  refine  the  angular  resolution  of  the  ray  fan  in  the 
model,  estimate  the  geo-acoustic  properties  of  the  bottom,  and  develop  transmission  loss 
estimates.  Transmission  loss  measurements  from  sono-buoy  data  were  used  as  an 
independent  metric  to  evaluate  the  model. 
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I.  INTRODUCTION 


A.  SOUTH  CHINA  SEA  OCEANOGRAPHY  AND  GEOMETRY 

The  South  China  Sea  (SCS)  is  a  marginal  sea  located  in  the  western  Pacific 
Ocean.  It  consists  of  a  deep  basin  surrounded  by  a  shallow  continental  shelf  to  the  north 
and  west,  and  tropical  islands  to  the  south.  The  deep  basin  portion  is  located  in  the 
eastern  central  portion  of  the  SCS  and  reaches  depths  in  excess  of  4,000  meters.  The 
SCS  is  connected  to  other  seas  and  the  north  Pacific  through  a  variety  of  straits,  the  most 
significant  of  which  is  the  Luzon  strait,  which  provides  the  great  majority  of  mass 
transport  into  the  SCS  (Ramp  et  al.  2004). 

The  SCS  is  home  to  significant  environmental  variability  on  scales  ranging  from 
the  super-tidal  to  seasonal,  which  consequently  induces  variability  in  the  SCS  acoustic 
propagation.  Some  of  the  more  interesting  features  include  variable  wind  forcing, 
intrusion  of  the  Kuroshio  current,  and  large  internal  waves.  The  SCS  wind  forcing 
features  summer  and  winter  monsoon  seasons,  as  well  as  tropical  cyclones.  The 
Kuroshio  current  is  the  strong  western  boundary  current  of  the  northwestern  Pacific 
ocean.  It  has  been  observed  to  enter  the  SCS  via  the  Luzon  strait,  and  its  subsequent 
behavior  is  widely  variable.  Possible  scenarios  include  bifurcation,  simply  exiting 
(meander),  or  forming  a  warm  core  ring  that  can  propagate  across  the  SCS.  (Ramp  et  al. 
2004)  Finally,  the  SCS  features  large  internal  waves  that  appear  to  be  generated  in  the 
Luzon  strait  via  the  abrupt  topography  change.  These  non-linear  waves  are  among  the 
largest  observed  and  can  impact  the  sound  speed  profile  and,  thus  the  acoustic 
propagation  in  the  basin.  One  of  the  major  goals  of  the  Windy  Island  Soliton  Experiment 
(WISE)  was  to  study  the  impact  of  the  internal  waves  on  acoustic  propagation. 

B.  WISE  EXPERIMENT  OVERVIEW 

The  Windy  Island  Soliton  Experiment  (WISE)  represents  a  collaborative  effort 
between  Taiwan  and  US  researchers  to  study  the  formation  and  evolution  of  the  trans¬ 
basin,  non-linear  internal  waves  in  the  Northeastern  South  China  Sea.  Moorings  were 
placed  on  the  continental  shelf,  the  deep  basin,  and  in  the  Luzon  Strait  for  a  period  of  a 
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year  commencing  in  April  2005.  In  addition  to  this  main  focus,  two  acoustic  propagation 
studies  were  conducted,  one  over  the  shallow  continental  shelf,  and  one  in  the  deep  basin 
in  order  to  study  the  physics  of  sound  propagation  in  the  SCS  environment,  and  the 
impact  of  the  large  internal  waves  on  sound  transmission.  In  particular,  the  focus  of  the 
basin  study  was  to  characterize  the  impact  of  the  trans-basin  internal  waves  on  long-range 
transmission  loss.  Physically,  this  was  accomplished  by  utilizing  two  deep-water 
moorings  that  transmitted  and  received  a  phase  modulated  400  Hz  signal  repeatedly  for  a 
period  of  a  year.  The  processed  signal  (the  output  of  a  correlation-matched  filter) 
corresponds  to  the  multipath  arrival  structure  of  a  pulse.  The  modeling  of  the  observed 
basic  arrival  structure  and  transmission  loss  in  a  “mean  ocean”  is  the  focus  of  this  thesis. 

C.  NAVAL  RELEVANCE 

The  Naval  relevance  of  this  thesis  takes  two  distinct  forms.  The  SCS  is  area  of 
significant  tactical  importance  to  the  Navy,  and  is  home  to  major  shipping  lanes  that 
account  for  a  significant  portion  of  world  commerce.  This  study  investigates  the 
observed  acoustic  arrival  structure  in  the  SCS  basin,  hopefully,  eventually  leading  to  an 
improved  understanding  of  acoustic  propagation  in  this  body  of  water.  Follow-on  studies 
will  focus  on  attempting  to  understand  the  acoustic  variability  as  a  function  of  trans-basin 
internal  waves,  which  will  provide  valuable  information  for  potential  operations.  In 
addition,  this  study  examines  the  geo-acoustic  properties  of  the  SCS  bottom,  valuable 
information  for  operational  planning  and  modeling. 

D.  APPROACH 

The  overall  objective  of  this  thesis  is  to  develop  and  refine  a  numerical  model  of 
the  400  Hz  signal  transmission,  in  particular  the  basic  arrival  structure  at  the  receiver  and 
the  transmission  loss  as  a  function  of  range  along  the  path.  Development  of  an  accurate 
model  of  the  basic  arrival  structure  and  transmission  loss  for  the  acoustic  path  used  in  the 
study  will  permit  future  modeling  study  on  the  impact  of  the  trans-basin  internal  waves 
on  sound  transmission  variability.  The  numerical  model  is  based  on  Hamiltonian  Ray 
tracing,  the  numerical  integration  of  Hamilton’s  equations  as  described  in  (Georges, 
Jones,  &  Riley,  1986)  and  (Chiu,  Semtner,  Ort,  Miller,  &  Ehret,  1994).  Known 
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environmental  data  includes  accurate  bathymetry  and  a  mean  sound  speed  profile 
determined  through  Conductivity-Temperature-Depth  (CTD)  casts  taken  during  the 
WISE  experiment.  There  are  three  primary  modeling  focuses  in  this  effort.  The  first  is 
an  empirical  detennination  of  the  necessary  angular  resolution  for  the  ray  fan. 
Equivalently,  over  a  given  angular  span  what  is  the  initial  ray  density  required  to 
accurately  model  the  observed  data.  The  second  focus  is  the  estimation  of  the  geo¬ 
acoustic  properties  of  the  sediment,  a  requirement  for  accurate  modeling  of  the 
amplitudes  of  the  multi-path  arrivals.  The  final  focus  is  the  development  of  a 
transmission  loss  curve  and  the  comparison  with  independent  measurements  as  an 
evaluator  of  the  model. 

E.  THESIS  OUTLINE 

The  remainder  of  this  thesis  consists  of  three  chapters,  Chapter  II  -  Methodology, 
Chapter  III  -  Modeling  and  Data  Comparison,  and  Chapter  IV  -  Results  and 
Conclusions.  An  appendix  containing  derivations  of  some  equations  is  also  included. 
Chapter  II  discusses  the  WISE  experiment  and  the  data  collected  in  more  detail  along 
with  Hamiltonian  Ray  tracing  and  the  Hamiltonian  Acoustic  Ray-tracing  Program  for  the 
Ocean  (HARPO)  model,  described  in  (Georges,  Jones,  &  Riley,  1986),  as  well  as  the 
approach  utilized  in  this  thesis  to  synthesize  the  arrival  structure.  Chapter  III  contains  the 
modeling  results  and  their  comparison  with  actual  data.  Chapter  IV  provides  a  summary 
of  significant  results. 
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II.  METHODOLOGY 


A.  WISE  TECHNICAL  PARAMETERS  AND  DATA 

1.  Basin  Experiment 

The  basin  acoustic  study  portion  of  the  WISE  experiment  consisted  primarily  of  a 
transmitter  and  receiver  moored  approximately  167  km  apart  with  the  source  and  receiver 
at  a  depth  of  approximately  840  m.  The  general  arrangement  is  depicted  on  Figure  1, 
where  the  locations  of  the  transmitter  and  receiver  are  designated  as  B2  and  B1 
respectively.  Additional  sono-bouy  measurements  were  taken  between  B2  and  B1  to 
measure  transmission  loss;  these  measurements  were  not  used  in  the  modeling  effort,  but 
rather  reserved  as  an  independent  evaluation  of  the  transmission  loss  predictions. 

The  sound  source  was  a  400  Hz  source  that  transmitted  a  5. 1 1  s,  5 1 1  digit  pseudo¬ 
random  signal  at  ~  1 83  dB  re  1  /iPa  at  1  m,  with  a  bandwidth  of  100  Hz.  The  source 
itself  consisted  of  a  six  m  array  of  four  hydrophones.  Signal  processing  gain  was 
achieved  by  pulse  compression  along  with  coherent  averaging  the  eight  m-sequence 
cycles  per  separate  reception.  The  final  signal  processing  gain  for  an  individual 
hydrophone  was: 

1 0  logI0  (8)  +  10  log10  (5 1 1)  =  36.1  dB 

The  desire  to  sample  the  internal  waves  in  the  basin,  along  with  hardware 
considerations  dictated  a  transmission  period  of  15  minutes.  The  moorings  were 
recovered;  the  batteries  were  serviced  and  redeployed  a  number  of  times  throughout  the 
year.  The  locations  and  transmitter/receiver  infonnation  used  are  from  the  final  such 
deployment.  Full  technical  details  of  the  entire  basin  experiment  are  available  in 
provided  in  (Miller  et  al.  2009) 
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Longitude  (Deg) 

Figure  1.  Locations  of  the  source  (B2)  and  receiver  (Bl)  in  the  SCS 
superimposed  on  the  bathymetry. 

2.  Data  Collected 

For  each  reception,  the  demodulated,  coherent  averaged  data  was  saved  for  each 
of  the  hydrophones.  The  data  that  was  utilized  in  this  study  consisted  primarily  of 
coherent  averaged  acoustic  arrival  structures  that  have  been  corrected  for  mooring  motion 
(Chiu  et  ah,  2008).  Some  examples  of  the  observed  arrival  structures  are  shown  in  Figure 
2.  It  is  clear  that  the  long  timeframe  of  this  experiment  admits  a  large  degree  of 
variability  into  the  individual  arrival  structures,  and  indeed  the  magnitude  of  the  first 
peak  varies  significantly,  as  can  be  seen  from  Figure  2.  While  the  intensity  of  the 
individual  peaks  in  the  arrival  varies  significantly,  the  basic  structure  itself,  which  is 
comprised  of  a  large  arrival  followed  by  a  number  of  smaller  arrivals,  is  fairly  consistent. 
Five  significant  peaks  were  consistently  observed  for  the  duration  of  the  WISE 
experiment,  as  is  evident  from  Figure  3.  Temporal  variability  is  clearly  present,  which  is 
not  the  subject  of  this  thesis.  Modeling  the  basic  arrival  structure  then  is  the  goal  of  this 
thesis,  along  with  a  reasonable  approximation  of  the  mean  intensity. 
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Figure  2.  Six  samples  of  acoustic  arrival  structures  observed  between  June  and 
October  2006. 
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Figure  3.  The  observed  arrival  structure  as  a  function  of  transmission  time  - 
displayed  in  dB  re  1  //  Pa  to  highlight  the  later  weak  arrivals,  (from 
Chiu  et  al.,  2008) 


B.  RAY  THEORY  AND  HAMILTONIAN  RAY  TRACING 

1.  General  Ray  Theory  and  the  Eikonal  Approximation 

In  general,  modeling  acoustic  propagation  through  the  ocean  necessitates  solving 
the  governing  wave  equation: 
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(2.1) 


1  82  p 

c(r )2  dt 2 


=  0 


Where  p  represents  pressure,  c  is  the  sound  speed,  t  denoted  time,  and  r  is  a  vector 
representing  position.  In  ray  acoustics  the  assumption  is  made  that  the  environmental 
change  and  thus  the  variation  in  sound  speed  is  gradual  over  a  small  number  of  acoustic 
wavelengths.  Ray  acoustics  involves  the  further  assumption  that  solutions  of  this 
equation  are  of  the  form: 


p(r,t)  =  A(r)  exp(i[2srft  -  k0fV(r)]) 


(2.2) 


2  nf 


Where  k0  =  — —  is  the  reference  wave  number,  and  W (r)  is  the  eikonal  and  represents 


the  phase  component  of  the  solution.  Since  solutions  of  constant  phase  represent  wave 
fronts,  and  rays  travel  in  a  direction  orthogonal  to  the  wave  front,  this  assumption  leads  to 
a  relation  that  governs  the  ray  paths.  The  ray  equations  are  given  below;  the  eikonal 
equation  yields  the  ray  equations  and  the  ray  paths  while  the  transport  equation  provides 


'c,  '  ' 


ray  amplitudes.  Here  the  ratio  of  sound  speed  to  a  reference  sound  speed  I  W 

\/c(r) 


is 


given  as  the  dimensionless  quantity  n ,  referred  to  as  the  index  of  refraction,  and  s  is  the 
distance  along  the  ray. 


Eikonal  Equation: 
Transport  Equation: 

Ray  Equations: 


2 


V  W  (r) 


n(r)2 


2VA(r)  •  V  W(r)  +  A(i)V2W(r)  =  0 

^r(n(r)^-)  =  Vn(r) 
as  as 


(2.3) 


A  full  derivation  of  these  equations  and  the  resultant  ray  equations  can  be  found  in  the 
Appendix  or  in  a  reference  such  as  (Ziomek  1995).  The  ray  approximation  is  also  known 
as  the  low  gradient  approximation  since  it  requires: 


V2A(r) 

K2ACr) 


1  <=>  Ix  f  » |  Vc 


(2.4) 
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For  a  transmission  frequency  of  400  Hz  the  quantity  2n  f  is  at  least  four  orders  of 
magnitude  larger  than|Vc|.  The  ray  paths  can  be  obtained  numerically  by  integrating 

Hamilton’s  equations,  which  are  equivalent  to  the  ray  equations  shown  above,  but 
derived  in  an  alternative  fashion  from  variational  principles. 


2.  Hamiltonian  Ray  Tracing 

Hamilton’s  equations  for  a  static  (no  current)  ocean  are  given  as: 

dr 

-j-  =  —VH(r,k) 
dr 

with  H(r,  k )  =  (In f  )2  -  c(r )2  \k2\ 


and  v  = - 

4  nf 


(2.5) 


Where  t  represents  travel  time  and  r  is  used  for  the  dependent  variable  in  Hamilton’s 
equations.  The  equivalence  is  based  on  a  variational  theme.  The  rays  travel 
perpendicular  to  the  wave  fronts  and  they  travel  along  the  path  of  shortest  time  subject  to 
the  constraint  of  the  index  of  refraction.  The  integration  of  Hamilton’s  equations  with  a 
suitably  defined  Hamiltonian  gives  just  such  variational  paths.  An  interesting  derivation 
that  uses  only  elementary  differential  operations  can  be  found  in  (Molcho  and  Censor 
1986).  An  overall  view  of  the  topic,  which  has  its  origins  in  geometrical  optics,  can  be 
found  in  (Combleet  1983)  where  Hamilton’s  approach  and  the  eikonal  approximation  are 
compared.  A  discussion  of  the  Hamiltonian  as  a  dispersion  relation  and  the  actual 
equations  integrated  is  given  in  the  Appendix. 


3.  Numerical  Integration 

This  approach  reduces  the  solution  of  the  free  wave  equation  on  some  section  of 

the  ocean  to  an  initial  value  problem  -  integrating  a  set  of  first  order  ODE’s.  The  initial 

conditions  specified  are  the  ray  launch  angles.  The  HARPO  code  in  FORTRAN 

originated  from  the  NOAA  Wave  Propagation  Laboratory  described  in  (Jones,  Riley  and 

Georges  1986)  It  integrates  Hamilton’s  equations  in  three  dimensions  in  polar  spherical 
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coordinates.  HARPO  utilizes  the  Adams-Moulton  predictor  corrector  for  the  numerical 
integration,  with  a  Runge-Kutta  starter.  The  starter  is  also  used  at  each  boundary 

interaction  and  turning  point.  The  algorithm  initially  sets  the  wave  number  k  equal  to 


the  value  expected  by  the  dispersion  relation,  that  is 


_(2 nfj_. 


c  00 


;  k  is  then  allowed  to 


vary  freely  as  Hamilton’s  equations  are  integrated.  The  error  is  then  evaluated  by  setting 


Error  = 


\k\ 


l*U 


-  - 1 ,  where  kdisp  is  the  wavenumber  expected  from  the  dispersion  relation. 


The  Runge-Kutta  (RK)  algorithm  is  utilized  to  provide  four  initial  points  to  start  the 
Adams-Moulton  (AM)  integration  step.  At  the  conclusion  of  each  AM  integration  the 
error  is  checked  against  a  user  provided  maximum  error.  If  the  error  has  been  exceeded, 
the  step  size  is  reduced  and  the  process  begins  again  with  the  RK  step.  This  approach 
has  the  advantage  of  some  numerical  efficiency  as  the  Adams-Moulton  method  is  an 
implicit,  multi-step  method  (of  0(h5))  that  varies  the  step  size  based  upon  the  error, 
which  is  related  to  the  gradient.  When  the  gradient  is  high,  the  step  size  will  be  low  and 
the  problem  will  be  more  computationally  complex.  When  the  gradient  is  low,  as  in  this 
modeling  effort,  the  step  size  can  be  large.  Full  details  of  the  particular  integrating 
algorithms  implemented  and  the  associated  errors  can  be  found  in  (Jones,  Riley  and 
Georges  1986). 


4.  Detecting  Bottom  and  Surface  Interactions 

HARPO  requires  bottom  data  that  is  continuous  and  smooth  enough  to  preclude 
very  irregular  shapes  where  diffraction  is  important  and  detection  of  bottom  interactions 
can  be  very  difficult.  HARPO  handles  bottom  and  surface  interactions  in  three  steps. 
First,  a  bottom  interaction  is  noted  at  some  integration  step  by  observing  that  a  surface 
(discontinuity  in  the  index  of  refraction)  has  been  crossed.  HARPO  then  iterates 
backwards  until  the  surface  is  reached,  and  finally  computes  a  reflection  direction  and 
then  begins  integrating  in  the  new  direction. 
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c. 


THE  HARPO  MODEL 


1.  Inputs  and  Technical  Data 

The  original  HARPO  model  described  in  (Georges,  Jones  and  Riley  1986) 
handled  only  canonical  oceans  and  boundary  conditions.  The  code  was  modified 
significantly  by  (Chiu  et  al.  1994)  to  permit  the  use  of  gridded  data  via  empirical 
orthogonal  eigenfunction  decomposition.  Bi-cubic  splines  are  fit  through  the  gridded 
data  to  create  continuous  representations  of  the  bottom  and  sound  speed  field.  This 
revision  also  created  MATLAB  code  that  interfaces  with  HARPO  for  ray  set  analysis, 
eigenray  search,  and  arrival  structure  synthesis. 

The  HARPO  model  is  fully  three  dimensional  and  can  support  three  dimensional 
sound  speed  data.  The  gridded  bathymetry  and  sound  speed  data,  along  with  the  source- 
receiver  geometry  are  the  required  data  inputs.  The  ray  launch  angles  are  the  initial 
conditions. 

2.  Ray  Set  Output 

The  output  of  the  HARPO  model  is  a  text  file  that  defines  the  ray  geometry.  It 
provides  the  ray  coordinates  in  spherical  polar  coordinates  referenced  to  the  mean  sea 
level.  This  of  course  can  be  imported  into  a  program  such  as  MATLAB  and  can  provide 
a  graphical  representation  of  the  ray  paths.  This  however  does  not  provide  an  arrival 
structure  or  any  method  to  determine  which  rays  interact  with  the  receiver.  Additionally, 
the  ray  set  provides  turning  points  and  interaction  with  the  surface  and  bottom.  This 
pennits  subsequent  analysis  to  consider  the  impact  of  bottom  interaction  on  the  energy 
transmitted  along  the  ray  path.  What  follows  describes  the  methodology  of  the  code 
developed  and  described  in  (Chiu,  Semtner,  Ort,  Miller,  &  Ehret,  1994). 

D.  RAYSET  ANALYSIS 

1.  Ray  Amplitude  Determination 

There  are  two  primary  factors  that  must  be  accounted  for  to  detennine  ray 
amplitude.  The  first  is  the  loss  of  energy  as  the  ray  travels  through  the  ocean.  In  ray 
acoustics  this  is  modeled  via  ray  tube  spreading,  the  concept  that  as  a  ray  tube  (bundle  of 
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rays)  travels  through  the  ocean  no  energy  is  lost  through  the  sides  of  the  tube  since  the 
rays  are  perpendicular  to  wave  fronts.  As  the  rays  travel  through  the  ocean,  the 
perpendicular  area  at  the  end  of  a  ray  tube  spreads  or  contracts,  resulting  in  a  different 
cross  sectional  area,  and  thus  a  varied  intensity.  This  quantity  is  inversely  proportional  to 
the  change  in  ray  range  with  respect  to  the  change  in  launch  angle  for  a  fixed  depth  and 
can  be  shown  to  be: 


P  =P  o 


2R^c_ 
r  cn 


cos  dn 


cos  0 


dz 


Rod0o 


(2.6) 


Where  p  is  pressure,  r  is  range  6{) ,  0  are  the  initial  launch  angle  and  the  ray  angle  at 


range  r  respectively.  R0  represents  the  reference  range  (1  m)  and  c0  the  sound  speed  at 

reference  range.  The  MATLAB  code  of  (Chiu,  Semtner,  Ort,  Miller,  &  Ehret,  1994)  is 
utilized  to  compute  this  derivative  numerically  via  finite  differences  using  adjacent  rays. 
The  code  utilizes  first  order  centered  differences  with  the  exception  of  the  rays  located  at 
the  edges  of  the  angular  span.  For  these  rays  first  order  forward  and  backward 
differences  are  used.  When  the  derivative  vanishes  (indicating  a  caustic)  or  the  ray  angle 

7T 

becomes  a  multiple  of—,  the  code  moves  outward  to  the  next  set  of  adjacent  rays.  This 


allows  the  code  to  track  through  caustics,  but  emphasizes  the  need  to  select  a  fine  enough 
angular  resolution  to  obtain  a  convergent  representation  of  the  ray  amplitude. 

The  second  factor  is  the  loss  in  ray  amplitude  due  to  an  interaction  with  the 
bottom  or  surface.  The  approach  utilized  here  is  that  of  (Clay  and  Medwin  1977),  where 
the  interaction  with  the  bottom  is  simulated  as  a  reflection  of  a  lossy  half  space.  A  full 
derivation  of  the  equations  given  below  is  located  in  the  aforementioned  reference.  The 
loss  is  accounted  for  by  considering  the  grazing  angle  of  the  ray  when  it  hits  the 
bottom  ,  the  water  column  sound  speed  c, ,  the  sediment  sound  speed  c2,  and  the  density 

of  the  water  column  and  bottom  px ,  p2 .  The  attenuation  of  compression  waves  in  the 
bottom  sediment  is  parameterized  by  letting  the  sound  speed  in  the  sediment  be  a 
complex  quantity  with: 
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(2.7) 


Im(c2)  = 


c2  af  1 
2nf  1000  201og10e 


This,  combined  with  the  assumption  of  a  plane  wave  solution  yields  the  following 
reflection  coefficient. 


p2  Re(c, )  sin  9X  +  i p2  Im(c2 )  sin  9X  -  pxcx  sin  02 
p2  Re(c2)sin  9X  +  ip2  Im(c2)sin  9X  +  pxcx  sin  92 


(2.8) 


Where  92  is  determined  from  Snell’s  law: 


92  =  arccosf 


Re(c2)  +  /'Im(c2) 


cos  9X  ] 


(2.9) 


The  MATLAB  code  of  (Chiu,  Semtner,  Ort,  Miller,  &  Ehret,  1994)  is  utilized  to  track  the 
ray  amplitude  through  each  interaction  with  the  bottom  and  surface. 


2.  Eigenray  Search  &  Parameters 

In  order  to  model  the  multipath  acoustic  arrival  structure  it  is  first  necessary  to 
identify  the  acoustic  paths  between  the  source  and  the  receiver,  the  so-called  eigenrays. 
The  MATLAB  code  of  Chiu,  Semtner,  Ort,  Miller,  &  Ehret,  1994)  utilized  to  search  the 
ray  set  and  determine  the  eigenray  parameters.  The  code  determines  the  depth  difference 
at  the  receiver  range  between  the  ray  and  the  receiver  for  each  ray  as  a  function  of  launch 
angle,  90 .  By  finding  the  zeroes  of  this  curve,  or  the  places  where  the  depth  difference  is 

zero,  the  precise  eigenray  parameters  can  be  obtained.  The  eigenray  determination  does 
not  involve  retracing  any  rays  -  for  the  purpose  of  plotting  eigenrays,  those  rays  nearest 
the  eigenrays  are  utilized.  A  sample  of  this  curve  is  shown  to  illustrate  the  process  in 
Figure.  4.  The  MATLAB  code  is  also  used  to  detennine  the  eigenray  parameters  such  as 
amplitude  an ,  phase  </>n ,  and  phase  front  arrival  time  tn . 
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depth  difference  vs  angle 


Figure  4.  A  curve  is  fit  through  the  depth  difference  vs.  launch  angle  for  each 
ray  -  the  zeroes  of  this  curve  represent  eigenrays  (shown  as  stars)  and 
permit  the  eigenray  parameters  to  be  precisely  determined. 

3.  Generating  Acoustic  Arrival  Structure 

Once  the  eigenray  parameters  have  been  determined,  an  acoustic  arrival  structure 
can  be  generated  via  the  following  equation.  A  model  pulse  is  created  to  match  the 
energy  and  shape  of  the  output  of  the  correlated  match  filter  used  in  signal  processing.  A 
comparison  of  these  two  pulses  is  shown  in  Figure  5.  Given  this  pulse,  and  the  eigenray 
parameters  identified  above,  the  resultant  acoustic  arrival  structure  can  be  expressed  as: 

r{t)=  Yj  a„4i-i„)expH(2^/oi„+^)]  (2.10) 

eigenrays 

Where  the  r,  s  indicated  the  complex  envelope  at  the  receiver  and  source 
respectively.  A  discussion  of  this  approximation  can  be  found  in  (Hager  2008). 
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Figure  5.  A  comparison  of  the  theoretical  pulse  used  to  synthesize  the  arrival 
structure  and  the  output  of  the  correlated  match  filter. 

4.  Modeling  Methodology 

The  next  chapter  details  the  application  of  this  methodology  to  the  specific  goals 
of  this  thesis.  First,  the  acoustic  arrival  structure  is  generated  at  increasingly  refined 
angular  resolutions  to  determine  what  resolution  is  necessary  to  obtain  a  convergent 
arrival  structure  while  minimizing  the  computational  intensity  of  the  problem.  Once  such 
a  satisfactory  ray  set  is  obtained  the  geo-acoustic  parameters  of  the  bottom  are  tuned  to 
match  the  structure  of  the  arrival  to  the  observed  data.  Finally,  transmission  loss 
predictions  are  generated  and  compared  to  independent  data  obtained  by  sonobuoys 
during  the  WISE  experiment. 
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III.  MODELING  AND  DATA  COMPARISON 


A.  BATHYMETRY  AND  SOUND  SPEED  DATA 

1.  Bathymetry 

The  two  sets  of  gridded  data  required  by  the  model  are  the  bathymetry  and  the 
sound  speed  profile.  The  bathymetry  utilized  was  extracted  from  a  high  resolution 
transect  from  the  Luzon  Strait  courtesy  of  (Liu  et  al.  1998).  The  input  data  is  referenced 
from  the  location  of  the  source  (B2)  and  reflects  the  deep  basin  nature  of  the  problem 
with  only  some  gentle  shoaling.  Some  smoothing  was  utilized  to  prevent  the  discrete 
nature  of  the  data  from  introducing  erratic  or  chaotic  behavior  in  the  ray  paths.  Figure  6 
shows  a  graphic  representation  of  the  bathymetry  with  the  source  and  receiver  annotated. 


Figure  6.  Gridded  data  input  to  HARPO,  the  sound  speed  profile  shown  in  the 
left  pane  is  a  mean  of  CTD  casts  taken  during  the  WISE  experiment 
and  the  bathymetry  shown  in  the  right  pane  is  that  between  the  source 
and  receiver. 

2.  Sound  Speed  Profile 

In  this  modeling  effort,  a  single  sound  speed  profile  (SSP)  was  utilized  to  create  a 
range  independent  sound  speed  field.  This  profile  represents  the  mean  of  the  CTD  casts 
taken  during  the  WISE  experiment.  The  location  of  the  source  and  receiver  in  the  SSP 
indicates  that  both  refracted-refracted  rays  and  bottom  interacting  rays  should  be 
expected.  Figure  7  shows  that  most  of  the  rays  interact  with  the  bottom  many  times, 
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making  the  bottom  composition  an  item  of  key  importance  in  modeling  the  arrival 
structure.  Some  smoothing  was  conducted  on  the  gridded  data,  as  well  as  extrapolating 
the  profile  in  the  gridded  data  set  to  ensure  that  it  extended  to  the  bottom  at  all  ranges  - 
ensuring  continuity  and  smoothness  of  the  bottom  and  sound  speed  representations. 

3.  Coordinate  System 

The  actual  latitude/longitude  coordinate  system  is  translated  to  a  simple  model 
coordinate  system  where  the  rays  can  be  shot  due  north.  The  ranges  and  curvature  of  the 
actual  coordinate  system  are  preserved.  The  HARPO  model  is  fully  capable  of  the 
complexities  of  a  three  dimensional,  range-dependent  problem,  the  focus  of  this  effort, 
however,  was  to  model  the  basic  arrival  structure  in  a  geometry  that  consisted  of  a  source 
and  receiver  in  line  with  each  other.  Therefore,  the  range  independent  sound  speed 
profile  and  the  bathymetry  were  duplicated  to  create  a  three-dimensional  box  where  the 
rays  will  behave  in  a  two  dimensional  fashion. 


Figure  7.  The  dependence  of  the  number  of  bottom  interactions  on  launch  angle 
-  showing  the  large  majority  of  rays  in  the  ray  fan  interact  with  the 
bottom  many  times. 
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B.  ANGULAR  RESOLUTION  STUDY 

1.  Angular  Aperture 

It  is  desirable  to  study  the  angular  resolution  to  ensure  that  no  possible  eigenpaths 
are  missed,  and  the  modeled  arrival  structure  converges  in  amplitude.  At  the  same  time, 
over  a  fixed  angular  aperture  an  increase  in  resolution  requires  a  larger  number  of  rays  to 
be  traced  which  in  turn  increases  the  computational  requirements.  It  is  desirable  then  to 
minimize  the  total  angular  aperture.  It  is  true  that  there  is  a  critical  angle  9o ,  beyond 

which  total  internal  reflection  occurs  -  but  this  is  a  function  of  the  bottom  and  since  a 
priori  we  cannot  be  certain  of  the  bottom,  we  must  select  an  angular  aperture  which  will 
contain  the  critical  angles  for  all  reasonable  bottom  choices.  The  geo-acoustic  properties 
of  the  bottom  and  the  characterization  of  the  bottom  in  this  effort  are  discussed  in  a  later 
section.  To  ensure  that  the  critical  angles  of  all  possible  bottoms  were  included  it  was 
necessary  to  utilize  an  aperture  of  at  least  [-20°,  20°  ] .  Figure  7  additionally  shows  that 
rays  outside  this  window  interact  with  the  bottom  many  times,  greater  than  ten,  and  thus 
could  not  contribute  to  the  arrival  structure  in  a  significant  way. 

2.  Angular  Resolution 

A  set  of  modeled  arrival  structures  resulting  from  a  progressive  refinement  in 
angular  resolution  is  shown  in  Figure  8.  The  angular  resolution  impacts  the  convergence 
of  the  arrival  structure  in  two  ways.  First,  it  is  necessary  to  ensure  that  all  contributing 
eigenpaths  are  captured.  Figure  8  shows  that  the  structure  is  fairly  constant  with  a 
resolution  as  course  as  0.05°,  indicating  that  at  this  resolution  all  contributing  eigenpaths 
have  been  captured.  The  second  concern  is  convergence  in  amplitude.  Since  the  ray  tube 
spreading  is  approximated  with  finite  differences,  it  is  necessary  to  utilize  an  angular 
resolution  fine  enough  to  obtain  a  sufficiently  convergent  representation  of  ray  amplitude. 
It  can  be  seen  from  Figure  8  that  the  amplitude  converges  at  the  finer  resolution  of 
O.Or/  0.005°.  Further  refinements  (0.0025°  /  0.001°)  confirmed  these  results.  The 
bottom  utilized  in  this  depiction  is  the  one  that  was  ultimately  selected  as  the  best  match: 
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a  silty-clay  bottom  whose  properties  are  discussed  below.  At  these  resolutions,  the  five 
arrivals  that  were  consistently  observed  during  the  WISE  effort  (Figure  3)  are  clearly 
observable. 


Time  (s) 


Figure  8.  Modeled  arrival  structures  corresponding  to  a  refinement  in  angular 
resolution;  the  structure  is  stable  at  0.05°,  but  convergence  in 
amplitude  requires  a  finer  resolution. 


3.  Eigenray  Composition 

The  eigenray  composition  of  each  of  the  five  observable  arrival  pulses  is  shown  in 
Figure  9.  The  first  arrival  occurs  at  approximately  112.6  seconds  and  is  comprised 


20 


entirely  of  refracted/refracted,  or  direct  rays.  The  second  arrival  occurs  at  1 12.7  seconds, 
and  is  comprised  of  eigenrays  that  interact  with  the  bottom  twice.  The  third  and  fourth 
arrivals  occur  at  112.7  and  112.9  seconds  respectively.  Here  the  number  of  bottom 
interactions  has  increased  to  three  and  four  bounces.  The  final  arrival  occurs  at  113.5 
seconds  and  consists  of  a  larger  number  of  eigenrays  that  interact  with  the  bottom  five 
times.  No  ray  that  interacts  with  the  bottom  more  that  this  (five  times)  contributes  to  a 
measurable  arrival. 
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Figure  9.  The  modeled  eigenray  composition  for  each  five  of  the  significant 
arrivals  that  were  consistently  observed  during  the  WISE  experiment. 
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c. 


BOTTOM  EVALUATION 


1.  Geoacoustic  Parameters 

The  ocean  bottom  acts  as  both  a  reflecting  and  scattering  surface  and 
characterizing  its  effect  on  ray  amplitude  is  a  complex  problem.  The  bottom  composition 
can  vary  from  rocky  to  muddy  surfaces.  Further,  even  a  specific  bottom  type  can  vary 
widely  in  terms  of  composition,  mean  grain  size,  porosity,  etc.  Often,  there  is  seasonal 
variability  caused  by  varying  sediment  deposition.  The  bottom  can  consist  of  multiple 
layers  of  different  material,  with  sound  penetrating  to  a  degree  based  on  its  frequency  and 
the  material  in  question.  The  problem  is  compounded  by  the  fact  that  in  deep  water, 
knowledge  of  the  bottom  composition  is  usually  limited  to  a  few  core  samples.  A  full 
discussion  of  the  acoustic  properties  of  the  sea  floor  can  be  found  in  (Hamilton  1980). 

The  HARPO  model  simulates  a  reflecting  ocean  bottom  and  tracks  the  bottom  and 
surface  interactions.  The  analysis  of  the  rayset  accounts  for  the  amplitude  loss  and  phase 
change  per  bottom  interaction  for  a  particular  eigenray  by  considering  the  grazing  angle 
of  the  ray,  the  water  column  sound  speed,  the  sediment  sound  speed,  the  density  of  the 
bottom,  and  the  attenuation  factor  a  dB  /  km  /  Hz ,  which  represents  the  attenuation  of 
compression  waves  by  the  bottom.  Some  more  sophisticated  models  differentiate 
between  this  attenuation  and  that  of  shear  waves,  but  this  refinement  is  not  considered  in 
this  effort.  In  general  a  is  determined  empirically.  Some  authors  (Urick  1983)  have 
noted  that  the  controlling  factor  in  the  acoustic  properties  seems  to  be  the  porosity  of  the 
sample,  but  (Richardson  and  Briggs  2004)  found  that  the  attenuation  can  depend  both  on 
the  intrinsic  attenuation  of  the  material,  and  that  caused  by  inhomogeneties  in  the 
substance,  and  thus  does  not  necessarily  regress  well  to  porosity.  Therefore  the  values  of 
a  utilized  in  this  modeling  effort  are  mainly  ones  that  have  been  determined  from  core 
samples  in  (Richardson  &  Briggs  2004)  or  postulated  by  (Schneck-Scott  2005). 

2.  Alternative  Cases 

Prior  to  the  conduct  of  the  WISE  experiment  modeling  was  conducted  as 
described  in  (Schneck-Scott  2005)  to  determine  the  resolvability  of  signals  in  the  SCS 
based  on  the  proposed  transmitter  and  receiver  arrangement.  Two  possible  bottom  cases 
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were  modeled  based  on  core  samples  that  indicated  a  bottom  that  ranged  from  sand-silt- 
clay  to  silty  clay.  These  two  possible  bottoms  were  intended  to  bind  the  problem  as  a 
best  and  worst  case.  In  (Richardson  &  Briggs  2004),  the  authors  survey  nearly  800  cores 
from  69  sites  around  the  world  in  a  compilation  of  physical  and  acoustic  properties. 
Cases  with  physical  properties  in  the  range  of  the  best  and  worst  case  bottoms  were 
selected  from  this  list  to  add  some  refinement  to  the  determination.  The  two  cases 
modeled  in  (Schneck-Scott  2005)  are  shown  in  Table  1.  The  alternative  cases  that  were 
examined  from  (Richardson  &  Briggs  2004)  are  shown  in  Table  2. 


Sediment  Type 

Sand-Silt-Clay 

Case  1 

Silty  Clay 

Case  2 

c  (m  /  s) 

1600 

1521 

P bottom  (kg /km3) 

1.7 

1.24 

a  (cIB  /  km  /  Hz) 

.313 

.486 

Table  1  Bottom  cases  modeled  in  (Schneck-Scott  2005)  in  preparation  for  the 
WISE  effort. 


Sediment  Type 

Clayey  Silt 

Case  A 

Sand-Silt 

Case  B 

Silty  clay 

CaseC 

Silty  Clay 
Case  D 

Silty  Clay 
Case  E 

c  (m  /  5) 

1545.5 

1614 

1507.2 

1480.4 

1501.7 

Pbonom  (kg  /km3) 

1.597 

1.820 

1.313 

1.506 

1.546 

a  (dB  /  km  /  Hz) 

0.579 

0.462 

0.285 

0.145 

0.166 

Table  2  Alternate  bottom  cases  drawn  from  (Richardson  &  Briggs  2004)  to 
provide  more  refinement  in  the  bottom  determination. 
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3. 


Bottom  Characterization 


A  comparison  of  the  modeled  arrival  structures  generated  by  the  bottoms  listed 
above  with  the  actual  arrival  structures  was  conducted.  It  was  possible  to  immediately 
reject  some  of  the  bottom  cases.  The  arrival  structures  generated  by  the  original  Case  1, 
and  the  alternative  cases  A,  and  B  were  extremely  inconsistent  with  the  observed  arrival 
structures,  containing  many  high  amplitude  peaks  spread  throughout  the  arrival  structure. 
The  structures  generated  by  Cases  2,  C,  D,  and  E  were  much  closer  to  the  “lossy”  bottom 
arrival  structures  that  were  observed.  Results  of  the  comparison  of  the  relative  amplitude 
of  the  peaks  indicated  that  Case  2,  the  original  worst  case  postulated  by  (Schneck-Scott, 
Sept,  2005)  was  the  closest  match  to  the  actual  bottom.  The  alternative  structures 
generated  by  Cases  C,  D,  and  E  were  in  fact  too  “lossy.”  The  structures  generated  by 
these  bottoms  are  shown  in  Figure  10,  where  it  can  be  seen  that  they  lack  the  smaller 
peaks  seen  the  actual  data.  Results  of  the  peak-to-peak  comparison  are  summarized  in 
Table  3,  for  completeness  all  cases  are  summarized.  The  bottom  then  is  determined  to  be 
silty-clay  bottom  characterized  as  in  Table  1,  the  original  Case  2. 
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Pressure  (Pa)  Pressure  (Pa)  Pressure  (Pa)  Pressure  (Pa) 


Alternate  Case  C 


Time  (s) 

Alternate  Case  D 


Time  (s) 

Alternate  Case  E 


Time  (s) 

Case  2  -  Closest  Match 


Figure  10  The  top  three  panes  show  the  reasonable  alternate  bottom  cases  - 
these  cases  lack  the  weak  late  arrivals  that  characterize  the  observed 
arrival  structure.  The  bottom  pane  shows  Case  2,  the  case  that  best 
matches  the  observed  data 
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Data  Set 

2nd  Peak 

3rd  Peak 

4th  Peak 

Actual 

-20  dB  Down 

-30-35  dB  Down 

-30-38  dB  Down 

Case  1 

-16  dB  Down 

-25  dB  Down 

-35  dB  Down 

Case  2 

-24  dB  Down 

-39  dB  Down 

-38  dB  Down 

Case  A 

-20  dB  Down 

-30  dB  Down 

-24  dB  Down 

Case  B 

-17  dB  Down 

-27  dB  Down 

-17  dB  Down 

CaseC 

-30  dB  Down 

-48  dB  Down 

-56  dB  Down 

Case  D 

-40  dB  Down 

-60  dB  Down 

-85  dB  Down 

Case  E 

-40  dB  Down 

-58  dB  Down 

-65  dB  Down 

Table  3  Comparisons  of  the  various  bottom  cases  to  observed  data  using  peak- 

to-peak  changes  on  dB  scale  referenced  from  the  first  peak;  Case  2 
provides  the  closest  match 


D.  MODELED  ACOUSTIC  ARRIVAL  STRUCTURE 

The  modeled  arrival  structure  generated  from  the  silty-clay  bottom  is  shown  at 
0.001°  resolution  in  Figure  8.  The  structure  compares  favorably  with  those  observed  in 
the  WISE  experiment,  with  the  peak  magnitude  well  within  the  observed  variability  of 
the  observations  made  during  the  WISE  experiment.  The  bottom  geometry  is  a 
dominating  factor  in  the  arrival  structure.  While  the  initial  peak  is  comprised  of 
refracted/refracted  (RR)  rays  -  the  subsequent  arrivals  are  comprised  of  eigenrays  that 
have  interacted  with  the  bottom  to  some  extent.  It  is  clear  that  the  vast  majority  of 
eigenrays  interact  with  the  bottom  numerous  times.  While  there  are  a  great  number  of 
eigenpaths,  only  a  select  few  of  them  actually  contribute  to  a  noticeable  arrival;  this  is  of 
course  due  to  the  lossy  nature  of  the  bottom. 
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E.  MODELED  TRANSMISSION  LOSS 
1.  Modeling  Transmission  Loss 

Transmission  loss  (TL)  is  an  ideal  tool  for  evaluation  of  acoustic  conditions  since, 
as  a  ratio  of  energies;  it  is  independent  of  source  strength.  Transmission  loss  is  defined 
as: 


TL  =  lOlog 


\p\(t)dt 

J  P2mM)dt 


(3-1) 


In  order  to  detennine  a  transmission  loss  curve  between  the  source  and  receiver,  a  model 
arrival  structure  was  generated  at  5  km  increments,  and  then  numerically  integrated  in 
accordance  with  equation  3.1.  This  analysis  was  conducted  at  a  model  depth  of  400  ft, 
since  that  was  the  depth  of  the  sonobuoy  measurements  taken. 


2.  Comparison  with  Sonobuoy  Data 

Figure  1 1  compares  the  TL  model  with  the  sonobuoy  data  taken  during  the  WISE 
experiment,  expressed  as  received  signal  intensity  level  (SL  -  TL) .  While  only  a  few 
sonobuoy  measurements  were  taken  in  the  region  modeled,  it  appears  that  they  agree  very 
well  with  the  model  prediction. 
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Recieved  Signal  Intensity  Level  (dES  re  (jPa) 


Received  Signal  Intensity  Level  -  Model  Prediction  vs  Sonobuoy  Measurements 


0  20  40  60  80  100  120  140  160 

Range  (km) 


Figure  11.  Comparison  of  the  modeled  received  level  with  the  sono-buoy  received 
levels  observed  during  the  WISE  experiment. 
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IV.  SUMMARY  AND  CONCLUSIONS 


The  objective  of  this  thesis  was  to  develop  and  refine  a  numerical  model  of  the 
400  Hz  transmission  utilized  in  the  SCS  WISE  experiment.  In  particular,  it  was  desired 
to  match  the  general  arrival  structure  observed  during  the  SCS  WISE  experiment.  This 
effort  utilized  a  Hamiltonian  ray-theory  approach  and  subsequent  analysis  of  the  ray  set 
in  order  to  model  the  acoustic  arrival  structure  generated  by  a  source  and  receiver  in  the 
SCS  basin.  A  mean  sound  speed  profile  developed  from  CTD  casts  was  utilized  to  create 
a  range-independent  sound  speed  field,  which  along  with  accurate  bathymetry  defined  the 
geometry  of  the  problem.  An  angular  resolution  study  was  conducted  to  establish 
convergence  of  the  arrival  structure  while  minimizing  computational  requirements.  The 
angular  aperture  for  the  rays  was  chosen  such  that  the  critical  angles  for  all  reasonable 
bottoms  fell  within  the  window  utilized.  This  study  established  that  an  angular  resolution 
of  0.005°  provided  sufficient  resolution  to  capture  all  contributing  eigenrays  and  provide 
a  sufficiently  convergent  estimate  of  the  amplitude  loss  due  to  ray  tube  spreading. 

Starting  with  a  general  notion  of  the  bottom  type  (silty-clay),  the  geo-acoustic 
properties  of  the  bottom  were  estimated  by  generating  acoustic  arrival  structures  for  a 
range  of  bottoms  and  comparing  them  with  the  WISE  observations.  With  both  the 
required  angular  resolution  and  the  bottom  parameters  determined,  a  transmission  loss 
model  could  be  estimated.  The  predicted  TL  agreed  very  well  with  the  sonobuoy 
measurements  that  were  conducted  during  the  WISE  experiment. 

Further  research  will  utilize  this  numerical  model  to  study  the  multi-scale 
variability  of  the  acoustic  arrival  structure,  and  the  impact  of  trans-basin  non-linear 
waves  on  acoustic  propagation.  This  future  research  will  attempt  to  understand  the 
physics  of  this  phenomena  and  the  impact  on  acoustics  in  the  SCS. 
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APPENDIX  -  DERIVATION  OF  THE  RAY  EQUATIONS 


A.l  The  Ray  Equations 

In  this  appendix,  a  derivation  of  the  ray  equations  is  provided.  The  approach 
taken  here  is  that  of  (Ziomek  1995).  A  full  discussion  of  the  ray  approach  can  be  found 
in  (Cornbleet  1983),  which  discusses  the  topic  through  the  lens  of  its  close  relative, 
geometric  optics.  Hamilton’s  approach  is  also  discussed  in  that  reference,  although  a 
simpler  approach  is  found  in  (Molcho  and  Censor  1986).  A  somewhat  qualitative 
description  of  ray  tracing  and  its  limitations  can  be  found  in  (Urick  1983). 


Beginning  with  the  free  wave  equation  and  substituting  a  solution  with  time 
hannonic  dependence,  we  obtain  the  Helmholtz  equation, 


where k0,  is  the  reference  wave  number  and  k0 


p{t,r)  =  A(r)eim 

=>  V2p(r)  +  k2n2(r)p(r)  -  0  (A.l) 

Where,  n  =  -^k- 
c(r ) 

The  further  assumption  that  the  solution  is  comprised  of  a  real  amplitude  part,  and  a  real 
phase  part , p(r)  =  A(r)e~,k°W(r) ,  where  W(r)  is  the  eikonal  discussed  in  II.B.l,  leads  to 
the  transport  equation  after  this  solution  is  substituted  into  the  Helmholtz  equation: 

[V2  +  k02n2(r)]p(r)  =  V2  (A(r)e~,k°w  <r))  +  k()2n2  (r)A(r)e~,k"ir<l)  =  0 

Evaluating  this  expression  with  vector  identities  and  equating  the  real  and  imaginary 
parts  with  zero  yields  the  system: 


V~A(r)  +  k02n2  (r)A(r)  -  k02A(r)  V  W(r) 


=  0 


2VA(r)  ■  V  W(r )  +  A(r)V2W(r )  =  0 


(A.2) 


The  second  of  these  equations  is  the  Transport  Equation  identified  in  Equation  2.3.  The 
first  of  these  equations  can  be  rewritten  as: 
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|2 


V2A(r)  |  1  |V^)| 
k2n2(r)A(r )  «2(r) 

or,  equivalently 


k2  VW(r)2 


V2A(r ) 

A(r) 


n2(r ) 


Here  we  assume  that  the  phase  is  changing  more  rapidly  than  the  amplitude,  or  rather  that 
the  change  in  the  amplitude  is  much  less  than  the  change  in  phase: 

k2VW(r)  »  — M  (A. 3) 

A(r) 

This  approximation  allows  the  first  equation  in  (A.2)  to  be  written  as: 

VW(r) 2  =  n2(r)  (A.4) 

This  equation  is  the  eikonal  equation  of  (2.3).  Recognizing  that  the  inequality  (A. 3)  is 
equivalent  to: 

V2A(r) 

K2A{r) 

yields  a  necessary  condition  of  ray  acoustics,  one  that  is  satisfied  in  this  modeling  effort 
as  discussed  in  Chapter  2. 

In  order  to  derive  the  ray  paths  from  the  eikonal  equation,  note  that  the  gradient  of  the 
eikonal  VW(r)  represents  the  direction  of  ray  propagation.  Then,  since  VW{r)  =  n(r), 
and  taking  the  derivative  of  both  sides  with  respect  to  s  ,  the  arc  length  along  the  ray,  we 
obtain: 

n  —  —  n  —dr 

-V(f(r)  =  V«(r)  =  -(n(r)-)  (A.5) 

ds  ds  ds 

This  coupled  set  of  ordinary  differential  equations  can  then  be  solved  to  obtain  the  ray 
paths. 


1  <=>  2  7i  f  » |Vc 
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A.2  The  Hamiltonian 


HARPO  calculates  ray  paths  by  integrating  Hamilton’s  equations.  In  general  the 
Hamiltonian,  H  is  defined  to  be  zero  along  ray  paths,  and  is  the  dispersion  relation  of  the 
acoustic  waves  set  equal  to  zero: 

H(r, k)  =  [a>-  hu{r)f  -  c 2  ( r)k 2  (A. 6) 

Here,  co  =  2nf ,  u(r)  is  the  current  field,  c(r)  is  the  sound  speed,  and  k  is  the  wave 
number.  In  this  modeling  effort,  a  static  (or  currentless)  ocean  is  simulated  and  (A. 6) 
reduces  to: 

H  =  or  -c2(r)k2  (A. 7) 

Physically,  the  quantity  C^->-  is  the  group  velocity  and  represents  the  speed  at  which 

dk 

energy  propagates  for  the  waves.  Since  HARPO  utilizes  polar-spherical  coordinates  the 

form  of  the  (A. 7)  utilized  is,  as  noted  in  (Georges,  Jones,  &  Riley,  1986): 

H  =  a 2  -  (kr2  +  kg2  +  k  2  )c2  (t,  r,  0,0)  =  0  (A. 8) 

In  order  to  obtain  the  system  of  ordinary  differential  equations,  the  partial  derivatives  of 
H  are  taken  with  respect  to  each  of  t,  r ,  0,  <j>,  kr  ,kg,k^.  This  yields  the  system  of  equations 

(Georges,  Jones,  &  Riley,  1986)  that  is  actually  numerically  integrated: 
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dH  ,  2  dc2 


dH 

dr 


- k 2 


dH 

~de 

dH 

d(j) 

dH 

dco 


=  -k 2 
=  -k 2 


=  2co 


dc -2_ 

dr 

dc2 

dd 

df_ 

d(j) 


dH  _ 
dkr 
dH  _ 

dkg 

dH  _ 

T  dH 
k>^^ 
dk 


-2  c2kr 
-lc2kg 
-2c\ 

=  -2c2 k2 
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